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g Abstract 

A highly efficient three-dimensional (3D) Lattice Boltzmann (LB) model for high speed com- 
pressible flows is proposed. This model is developed from the original one by Kataoka and Tsuta- 
hara[Phys. Rev. E 69, 056702 (2004)]. The convection term is discretized by the Non-oscillatory, 
\ containing No free parameters and Dissipative (NND) scheme, which effectively damps oscillations 

at discontinuities. To be more consistent with the kinetic theory of viscosity and to further improve 
the numerical stability, an additional dissipation term is introduced. Model parameters are chosen 
in such a way that the von Neumann stability criterion is satisfied. The new model is validated 
\ by well-known benchmarks, (i) Riemann problems, including the problem with Lax shock tube 

t^J- ■ and a newly designed shock tube problem with high Mach number; (ii) reaction of shock wave on 

d 



O 



droplet or bubble. Good agreements are obtained between LB results and exact ones or previ- 
ously reported solutions. The model is capable of simulating flows from subsonic to supersonic and 
capturing jumps resulted from shock waves. 
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I. INTRODUCTION 



Lattice Boltzmann (LB) method has been becoming a powerful and efficient tool to 



simulate fluid flows in many areas fl], rangin 
dynamics (46], flows through porous media [ 



rom multiphase flows 



, y| , magnetohydro- 



and thermal fluid dynamics [9] . However, 



most models so far work only for incompressible fluids. Many attempts have been made in 



constructing LB models for the compressible Euler equations. Hu et al. [10] proposed a 
13-discrete-velocity model based on the triangular lattice. In this model, particles at each 
node are classified into three kinds. They are on the energy levels e^, e#, and e^, where 
£a > £b > 0, the energy level en is higher than and is for the rest particle. Similar to Hu's 



model, Yan and co-workers 



ll| presented a 17-discrete- velocity model with three-speed- 



three-energy level on a square lattice. Both models are two-dimensional (2D) and belong to 
the standard LB model. In the standard LB model, particle velocities are restricted to those 
exactly linking the lattice nodes in unit time. Besides the standard LB, Finite Difference 
(FD) LB is attracting more attention with time. With the FD LB model we do not need 
consider that constraint, we can choose particle velocities independently from the lattice 
configuration. 

n 

Shi et al. [12j formulated a FD LB scheme based on a two-dimensional 9-velocity model. 
This model allows particles to possess both kinetic and thermal energies. Kataoka and 
Tsutahara [l3| presented a LB model series for the compressible Euler equations, where 
5, 9 and 15 discrete velocities are used for the one- , two- and three-dimensional cases, 
respectively. However, all these models work only for subsonic flow. The low-Mach number 
constraint is generally related to the numerical stability problem. The latter has been partly 
addressed by a few potential solutions, for example, the entropic method [ijj], flux limiters 

and multiple-relaxation-time LB approach j^ . 



dissipation techniques 



16 



Watari and Tsutahara proposed a three-dimensional FD LB model for Euler equations, 
where numerical simulations are successfully performed up to Mach number 1.7 2lj. But 
the number of discrete velocities in that model is up to 73, which is quite expensive from 
the view of computational side. Recently, a three-dimensional compressible FD LB model 
without free parameters was proposed [22J, where 25 discrete velocities are used. With 
this model the momentum equations at the Navier-Stokes level and energy equation at the 
Euler level can be recovered. The maximum Mach number is 2.9 in simulations. Pan, et 
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al. 



16| developed the 2D model by Kataoka and Tsutahara [13| by introducing reasonable 



dissipation term so that the model works for supersonic flows. Flows with Mach number 
higher than 30 are successfully simulated with the model. 

In this paper we formulate a three-dimensional FD LB model for high speed compressible 
flows, based on Kataoka's 15-velocity model and reasonable dissipation technique. The 
following part of the paper is planned as follows. Section 2 presents the discrete velocity 
model used in this work. Section 3 describes briefly the FD scheme and performs the von 
Neumann stability analysis. Simulation results are presented and analyzed in Section 4. 
Section 5 makes the conclusion. 



II. 3D DISCRETE VELOCITY MODEL BY KATAOKA AND TSUTAHARA 



The evolution of the distribution function f\ is governed by the following equation 



23|: 



where Vi a is the a component of velocity i>j, % = 1,...,N, N is the number of discrete 

velocities, index a = 1, 2, 3 corresponding to x, y, and z, respectively. The Einstein's 

convention for sums is used. The variable t is time, x a is the spatial coordinate, f^ q is 

the local-equilibrium distribution function, and r represents the relaxation time. At the 

continuous limit, the above formulation is required to recover the following Euler equations: 

dp_ djpug) = Q 
dt dx a 

d(pug) + d(pu a up) + &P_ = Q , 2 x 
dt dxp dx a 

dp(bRT + u 2 a ) | dpu a (bRT + ul) + 2Pu a _ q 
dt dxp 

where p, u a , T, P are, respectively, the density, the flow velocity in the x a direction, the 

temperature, and the pressure of gas. R is the specific gas constant and b is a constant 

relating to the specific-heat ratio 7, b = 2/(7 — 1). The 3D discrete velocity model proposed 

by Kataoka and Tsutahara (see Fig. 1) can be expressed as: 



{Vil,Vi2,V i3 ) 



(0, 0, 0) for i = 1, 

ci (±1, 0, 0) , ci (0, ±1,0), ci (0, 0, ±1) for % = 2, 3, • • • , 7, 

^(±1,±1,±1) for 1 = 8,9,- •• ,15, 
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FIG. 1: Distribution of Vj for the proposed discrete velocity model. 



^ = \ (3) 
0, fori = 2,3,--- ,15, 

where Ci,c 2 , and 7y are given nonzero constants. In this model, the local-equilibrium distri- 
bution function f^ q satisfies the following relations: 

TV 

P = £/f, (4a) 
1=1 

TV 

P«a = J] (4b) 

1=1 

TV 



p(6i?T + ^) = ^/r(^ + ^), (4c) 
1=1 

TV 

P<J aj a + pM a M/3 = ^ fi^^ifi, ( 4d ) 



i=l 
TV 



p[(b + 2)iZT + u%u a = J2 ft>l + Vi)vic (4e) 

i=i 

The local-equilibrium distribution function is defined as follows: 

/r = P(A + BiV ia U a + DiUaViallpVip), i = l,2,..., 15, (5) 
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where 



A; 



B: 



1 

6(c?-eg) 
1 

^ 8(3 -eft 

0. 



'/d 



((6-3)4 + 3)r + > 2 



-c?+((6-3)| + 3)T + ^ 



2c|" 



C L,2 



- C 2+(6+2)T+«2 
3[- C 2 + (6+2)T+n2] 



1 

2,3,- 
= 8,9, 



15 



Di 



16(4 



= 1 

2,3, 



8,9, 



= 2,3, 
8,9, 



15 



,7 
,15 



(6) 



III. FD SCHEME AND VON NEUMANN STABILITY ANALYSIS 

In the original LB model the finite difference scheme with the first-order forward 
in time and the second-order upwind in space is used for the numerical computation. This 
model has been validated via the Riemann problem in subsonic flows and encounters in- 
stability problems in supersonic flows. In order to improve the stability, we adopt the 
Non-oscillatory, containing No free parameters and Dissipative (NND) scheme for space dis- 
cretization. To be more consistent with the kinetic theory of viscosity and to further improve 
the numerical stability, an additional dissipation term is introduced. 

In the NND scheme, the spacial derivative is calculated using the following formula: 

d(v ia fi) 1 



where / represents node index in x or y direction. is the numerical flux at the interface 

of (xi + 4p, y) or (x, yi + ^) , and defined as: 



Ku 



L I fR 



(8) 



where 



f L T ^i = f, + r + ~ minmod ( A f+ x , A f+ , 
f R r , i = fTru-1 — - min mod ( A f ~ , , A f ~ 

' \Via\) fi,I, 



Af ± , = f 



>i,I+l 



minmod(X,F) = -min(|X| , \Y\) [Sign (X) + Sign (Y)] . 



(9) 



The NND scheme itself contains a forth-order dissipation term with a negative coefficient 
which reduces the oscillations, but it is not enough to highly improve the stability, which 
means an additional dissipation term is needed for a practical LB simulation. In order to 
further improve the stability, and enhance its applicability for high Mach flows, we introduce 
artificial viscosity into the LB equation: 

t + ^ = ~<'<-^> + *£g. (io) 

a=l " 

where 

ciAx, i = l 

K = < a Ax/10, i = 2,3, ••• ,7 . 

0, z = 8,9, •-• ,15 

The second-order derivative can be calculated by the central difference scheme. 

In the following we do the von Neumann stability analysis of the improved LB model. In 

the stability analysis, we write the solution of FD LB equation in Fourier series form. If all 

the eigenvalues of the coefficient matrix are less than 1, the algorithm is stable. 

Distribution function is split into two parts: fi(x a ,t) = ff + Afi(x a ,t), where is the 

global equilibrium distribution function. It is a constant which does not change with time 

or space. Putting this equation into Eq. (fTDl) we obtain: 

Afj(x Q ,t + At) - Afj(x a ,t) djj_ _ _1 _ eq . (Pfi 

At + Vta dx a ~ r [ h h ' + Al dxl ' { ' 

the solution can be written as 

Afi(x a ,t) = F*exp(ik a x a ), (12) 

where F\ is an amplitude of sine wave at lattice point x a and time t, k a is the wave number. 
From the Eq.tfTTl) and Eq.( fl2l) we can get F- +At = GijFj. Coefficient matrix Gij describes 
the growth rate of amplitude F\ in each time step At. The von Neumann stability condition 
is max|o;| < 1, where u denotes the eigenvalue of coefficient matrix. Coefficient matrix 
of NND scheme can be expressed as follows, 

/ At V ia At\ x Atdf? ( e ik a A Xa _ 2 + e -ik a A Xa ) 

G « = ~ T ~ -Ax^V S « + ~M 1 i^J 2 ^ 



(1 - 5) (1 - e~ ikaAxa ) , if v ia > 0; 
(l - ^) (e ikaAx * - 1) , if v ia < 0. 



1.0010- 
1.0005- 
1.0000- 

s 

£ 0.9995- 
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0.9990 - 
0.9985 - 
0.9980 - 



+++++ 



+ - + ;° oo % 
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kdx 



FIG. 2: Stability analysis of several finite difference schemes. 

dfj dp <).(, or <).(, du a oj) ■ ° 2' 

There are some numerical results of von Neumann stability analysis by Mathematica. 
Abscissa is kdx, and ordinate is Iculm^ that is the biggest eigenvalue of coefficient matrix 
Gij. 

Figure 2 shows the stability analysis of several finite difference schemes. The macroscopic 
variables are set as (p, u 1: u 2 , u 3 , T) = (1.0, 4.0, 0.0, 0.0, 1.0), the other model parameters are: 
( Cl ,c 2 ,Vo) = (4.0,12.0,4.0), dx = dy = dz = 4 x 10~ 3 , dt = r = 10~ 5 , b = 5, a = p = 0. 
In this test, the NND scheme shows better stability than the others. Figure 3 shows the 
effect of dissipation term. The variables are set as (p, u±, w 2 , w 3 , T) = (1.0, 20.0, 0.0, 0.0, 1.0), 
(ci,C2,?7o) = (20.0,60.0,20.0), and the others are consistent with the Figure 2. In the two 
cases of Figure 3, operation with dissipation term is more stable (max \u\ < 1). 

Figure 4 shows the influence of parameters c±, c 2 , i] on the stability in the absence of 
dissipation term. The macroscopic variables and the other model parameters are consistent 
with those of Figure 2. Figure 5 shows the stability effect of the three parameters, when 
there is a dissipation term. The macroscopic variables and the other model parameters are 
consistent with those of Figure 3. In Figure 4 constants Ci, c 2 and f] affect the stability 
heavily. In Figure 5 the LB is stable for all tested values of c 2 and i] . Based on these tests, 
we suggest that C\ can be set a value close to the maximum of flow velocity, c 2 can be chosen 
about 3 times of the value of c±, and r] can be set to be about 1 ~ 2 times of the value of 
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< 2- 



(13) 
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FIG. 3: Effect of dissipation term on numerical stability. 




FIG. 4: Influence of parameters c±, C2, f]o on stability in the absence of artificial viscosity. 
IV. NUMERICAL SIMULATION AND ANALYSIS 

In this section we study the following questions using the proposed LB model: one- 
dimensional Riemann problems, and reaction of shock wave on a droplet or bubble. 
(1) One-dimensional Riemann problems 

Here, we study two one-dimensional Riemann problems, including the problem with Lax 
shock tube and a newly designed shock tube problem with high Mach number. Subscripts 
"L" and "R" indicate the left and right macroscopic variables of discontinuity. 

(a) Lax shock tube problem 
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FIG. 5: Effect of c±, c 2 , rjo under the condition with artificial viscosity term. 

The initial condition of the problem can be defined: 

(p, Ul ,u 2 ,u 3 ,T)\ L = (0.445,0.698,0.0,0.0,7.928), 
(p, Ul ,u 2 ,u 3 ,T)\ R = (0.5, 0.0, 0.0, 0.0, 1.142). 

Figure 6 shows the comparison of the NND scheme and the second-order upwind scheme 
without the dissipation term at t — 0.1. Circles are for the NND scheme simulation results, 
squares correspond with the second-order upwind scheme, and solid lines are for exact so- 
lutions. The parameters are (ci,c 2 ,r]o) = (2.0,6.0,2.0), 7 = 1.4, dx = dy = dz = 0.003, 
dt = t = 10 -5 . Compared with the simulation results of second-order upwind scheme, the 
oscillations at the discontinuity are weaker in the NND simulation, 
(b) High Mach number shock tube problem 

In order to test the Mach number of the new model, we construct a new shock tube 
problem with high Mach number, and the initial condition is 

(p,u 1 ,u 2 ,u 3 ,T)\ L = (100.0,10.0,0.0,0.0,0.714286), 

(15) 

(p,u u u 2 ,u 3 ,T)\ R = (150.0,0.0,0.0,0.0,50.0). 

Figure 7 shows a comparison of the numerical results and exact solutions at t — 0.25, 
where (c 1 ,c 2 ,r]o) = (8.0,24.0,8.0), 7 = 1.4, dx = dy = dz = 0.01, dt = r = 10~ 5 . The 
Mach number of the left side is 10 (Ma = u/^jT = 10/y/lA x 0.714286), and the right is 
(Ma = -u/VtT = 0). Successful simulation of this test shows the proposed model is still 
likely to have a high stability when the Mach number is large enough. 

9 




FIG. 6: Numerical results and exact solutions for Lax shock tube at t = 0.1. 
(2) Reaction of shock wave on 3D bubble problem 

The proposed model is used to simulate interaction of a planar shock wave with a bubble 
or droplet. The shock wave is moving from the right to the left. Initial conditions are (a) 



I 



(1,0,0,0,1), pre — shock, 

(p,ui,u 2 ,u 3 ,p) \ Xt y t0 = { (2.66667,-1.47902,0,0,3.94406), post - shock, (16) 

(0.1358,0,0,0,1), bubble, 



and (b) 



(p,u 1 ,u 2 ,u 3 ,p) \ x>vfl = 



(1,0,0,0,1), pre — shock, 

(2.66667,-1.47902,0,0,3.94406), post - shock, (17) 
(4.1538,0,0,0,1), bubble. 



The corresponding shock wave Mach number is 2.0, (Ma = (D — u) / V7T = (2.36643 — 
0)/Vl-4 x 1, where D = 2.36643 is the wavefront velocity). 

The domain of computation is (0 : 301,0 : 81,0 : 81). Initially, the bubble or droplet is 
at the position (200,40,40). In the simulations, the right side adopts the values of the initial 
post-shock flow, the extrapolation technique is applied at the left boundary, and reflection 
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FIG. 7: The numerical and exact solutions for high Mach number shock tube at t = 0.25. 



conditions are imposed on the other four surfaces. Specifically, at the right side, 

p(NX + 1, iy, iz) = p{NX, iy, iz) = 2.66667, 
T(NX + 1, iy, iz) = T(NX, iy, iz) = 1.6875, 
u^NX + l,iy,iz) = m(NX,iy,iz) = -1.47902, 

u 2 (NX + 1, iy, iz) = u 2 (NX, iy, iz) = 0, 

Us(NX + 1, iy, iz) = Us(NX, iy, iz) = 0, 

where ix (or iy, iz) is the index of lattice node in the x- (or y-, z-) direction, and ix = 0, 
1, • • • , NX + 1 ( iy = 0, 1, • • • , NY + 1; iz = 0, 1, • • • , NZ + 1). At the left side 
p{l,iy,iz) = 2p(2,iy,iz) — p(3,iy,iz), p(0,iy,iz) = 2p(l,iy,iz) — p(2,iy,iz), temperature 
and velocity components have the same form. Finally we take the upper surface as an 
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example to describe the reflection conditions. 



p{ix,NY + l,iz) 
T(ix, NY + l,iz) 
m(ix,NY + l,iz) 
U2{ix, NY + 1, iz) 
u^{ix, NY + 1, iz) 



= p{ix, NY — 1, iz), 

= T(ix, NY - l,iz), 

= Ui(ix, NY — 1, iz), 

= — u<i{ix, NY — 1, iz), 

= u^iix, NY — l,iz), 



p(ix, NY, iz) 
T(ix, NY, iz) 
Ui{ix, NY, iz) 
u 2 (ix, NY, iz) 
Us(ix, NY, iz) 



p{ix, NY — 1, iz), 
T(ix, NY — 1, iz), 
Ui{ix, NY — 1, iz), 
0. 

Us(ix, NY — l,iz). 



Parameters are as follows: (ci,c 2 ,?7o) — (2.0,6.0,4.0), 7 = 1.4, dx = dy = dz = 0.001, 
dt = t = 10~ 5 . Figure 8 and Figure 9 show the density iso-surfaces of bubble or droplet, 
where Figure 8 is for the process with initial condition (Tl6|) . and Figure 9 is for condition 
f fT7|) . Figure 10 shows the density contours on section z = 40, where (a) and (b) correspond 
to the processes of Figure 8 and Figure 9, respectively. The simulation results are accordant 
with those by other numerical methods 24j, |25[ and experiment 26 1. 



V. CONCLUSION 



We proposed a highly efficient 3D LB model for high-speed compressible flows. The 
convection term in Boltzmann equation is solved with the finite difference NND method, 
additional dissipation term is introduced to match the more realistic kinetic viscosity and 
to be more stable in numerical simulations. Model parameters are controlled in such a way 
that the von Neumann stability criterion is satisfied. The model can be used to simulate 
flows from subsonic to supersonic flows, especially supersonic flows with shock waves. 
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y 

FIG. 8: Density iso-surfaces of a low density bubble at t = 0.0, 0.1, respectively. 
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I. INTRODUCTION 



Lattice Boltzmann (LB) method has been becoming a powerful and efficient tool to 



simulate fluid flows in many areas fl], rangin 
dynamics (46], flows through porous media [ 



rom multiphase flows 



, y| , magnetohydro- 



and thermal fluid dynamics [9] . However, 



most models so far work only for incompressible fluids. Many attempts have been made in 



constructing LB models for the compressible Euler equations. Hu et al. [10] proposed a 
13-discrete-velocity model based on the triangular lattice. In this model, particles at each 
node are classified into three kinds. They are on the energy levels e^, e#, and e^, where 
£a > £b > 0, the energy level en is higher than and is for the rest particle. Similar to Hu's 



model, Yan and co-workers 



ll| presented a 17-discrete- velocity model with three-speed- 



three-energy level on a square lattice. Both models are two-dimensional (2D) and belong to 
the standard LB model. In the standard LB model, particle velocities are restricted to those 
exactly linking the lattice nodes in unit time. Besides the standard LB, Finite Difference 
(FD) LB is attracting more attention with time. With the FD LB model we do not need 
consider that constraint, we can choose particle velocities independently from the lattice 
configuration. 

n 

Shi et al. [12j formulated a FD LB scheme based on a two-dimensional 9-velocity model. 
This model allows particles to possess both kinetic and thermal energies. Kataoka and 
Tsutahara [l3| presented a LB model series for the compressible Euler equations, where 
5, 9 and 15 discrete velocities are used for the one- , two- and three-dimensional cases, 
respectively. However, all these models work only for subsonic flow. The low-Mach number 
constraint is generally related to the numerical stability problem. The latter has been partly 
addressed by a few potential solutions, for example, the entropic method [ijj], flux limiters 

and multiple-relaxation-time LB approach j^ . 



dissipation techniques 
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Watari and Tsutahara proposed a three-dimensional FD LB model for Euler equations, 
where numerical simulations are successfully performed up to Mach number 1.7 2lj. But 
the number of discrete velocities in that model is up to 73, which is quite expensive from 
the view of computational side. Recently, a three-dimensional compressible FD LB model 
without free parameters was proposed [22J, where 25 discrete velocities are used. With 
this model the momentum equations at the Navier-Stokes level and energy equation at the 
Euler level can be recovered. The maximum Mach number is 2.9 in simulations. Pan, et 



2 



al. 



16| developed the 2D model by Kataoka and Tsutahara [13| by introducing reasonable 



dissipation term so that the model works for supersonic flows. Flows with Mach number 
higher than 30 are successfully simulated with the model. 

In this paper we formulate a three-dimensional FD LB model for high speed compressible 
flows, based on Kataoka's 15-velocity model and reasonable dissipation technique. The 
following part of the paper is planned as follows. Section 2 presents the discrete velocity 
model used in this work. Section 3 describes briefly the FD scheme and performs the von 
Neumann stability analysis. Simulation results are presented and analyzed in Section 4. 
Section 5 makes the conclusion. 



II. 3D DISCRETE VELOCITY MODEL BY KATAOKA AND TSUTAHARA 



The evolution of the distribution function f\ is governed by the following equation 



23|: 



where Vi a is the a component of velocity i>j, % = 1,...,N, N is the number of discrete 

velocities, index a = 1, 2, 3 corresponding to x, y, and z, respectively. The Einstein's 

convention for sums is used. The variable t is time, x a is the spatial coordinate, f^ q is 

the local-equilibrium distribution function, and r represents the relaxation time. At the 

continuous limit, the above formulation is required to recover the following Euler equations: 

dp_ djpug) = Q 
dt dx a 

d(pug) + d(pu a up) + &P_ = Q , 2 x 
dt dxp dx a 

dp(bRT + u 2 a ) | dpu a (bRT + ul) + 2Pu a _ q 
dt dxp 

where p, u a , T, P are, respectively, the density, the flow velocity in the x a direction, the 

temperature, and the pressure of gas. R is the specific gas constant and b is a constant 

relating to the specific-heat ratio 7, b = 2/(7 — 1). The 3D discrete velocity model proposed 

by Kataoka and Tsutahara (see Fig. 1) can be expressed as: 



{Vil,Vi2,V i3 ) 



(0, 0, 0) for i = 1, 

ci (±1, 0, 0) , ci (0, ±1,0), ci (0, 0, ±1) for % = 2, 3, • • • , 7, 

^(±1,±1,±1) for 1 = 8,9,- •• ,15, 

3 



FIG. 1: Distribution of Vj for the proposed discrete velocity model. 



^ = \ (3) 
0, fori = 2,3,--- ,15, 

where Ci,c 2 , and 7y are given nonzero constants. In this model, the local-equilibrium distri- 
bution function f^ q satisfies the following relations: 

TV 

P = £/f, (4a) 
1=1 

TV 

P«a = J] (4b) 

1=1 

TV 



p(6i?T + ^) = ^/r(^ + ^), (4c) 
1=1 

TV 

P<J aj a + pM a M/3 = ^ fi^^ifi, ( 4d ) 



i=l 
TV 



p[(b + 2)iZT + u%u a = J2 ft>l + Vi)vic (4e) 

i=i 

The local-equilibrium distribution function is defined as follows: 

/r = P(A + BiV ia U a + DiUaViallpVip), i = l,2,..., 15, (5) 



4 



where 



A; 



B: 



1 

6(c?-eg) 
1 

^ 8(3 -eft 

0. 



'/d 



((6-3)4 + 3)r + > 2 



-c?+((6-3)| + 3)T + ^ 



2c|" 



C L,2 



- C 2+(6+2)T+«2 
3[- C 2 + (6+2)T+n2] 



1 

2,3,- 
= 8,9, 



15 



Di 



16(4 



= 1 

2,3, 



8,9, 



= 2,3, 
8,9, 



15 



,7 
,15 



(6) 



III. FD SCHEME AND VON NEUMANN STABILITY ANALYSIS 

In the original LB model the finite difference scheme with the first-order forward 
in time and the second-order upwind in space is used for the numerical computation. This 
model has been validated via the Riemann problem in subsonic flows and encounters in- 
stability problems in supersonic flows. In order to improve the stability, we adopt the 
Non-oscillatory, containing No free parameters and Dissipative (NND) scheme for space dis- 
cretization. To be more consistent with the kinetic theory of viscosity and to further improve 
the numerical stability, an additional dissipation term is introduced. 

In the NND scheme, the spacial derivative is calculated using the following formula: 

d(v ia fi) 1 



where / represents node index in x or y direction. is the numerical flux at the interface 

of (xi + 4p, y) or (x, yi + ^) , and defined as: 



Ku 



L I fR 



(8) 



where 



f L T ^i = f, + r + ~ minmod ( A f+ x , A f+ , 
f R r , i = fTru-1 — - min mod ( A f ~ , , A f ~ 

' \Via\) fi,I, 



Af ± , = f 



>i,I+l 



minmod(X,F) = -min(|X| , \Y\) [Sign (X) + Sign (Y)] . 



(9) 



The NND scheme itself contains a forth-order dissipation term with a negative coefficient 
which reduces the oscillations, but it is not enough to highly improve the stability, which 
means an additional dissipation term is needed for a practical LB simulation. In order to 
further improve the stability, and enhance its applicability for high Mach flows, we introduce 
artificial viscosity into the LB equation: 

t + ^ = ~<'<-^> + *£g. (io) 

a=l " 

where 

ciAx, i = l 

K = < a Ax/10, i = 2,3, ••• ,7 . 

0, z = 8,9, •-• ,15 

The second-order derivative can be calculated by the central difference scheme. 

In the following we do the von Neumann stability analysis of the improved LB model. In 

the stability analysis, we write the solution of FD LB equation in Fourier series form. If all 

the eigenvalues of the coefficient matrix are less than 1, the algorithm is stable. 

Distribution function is split into two parts: fi(x a ,t) = ff + Afi(x a ,t), where is the 

global equilibrium distribution function. It is a constant which does not change with time 

or space. Putting this equation into Eq. (fTDl) we obtain: 

Afj(x Q ,t + At) - Afj(x a ,t) djj_ _ _1 _ eq . (Pfi 

At + Vta dx a ~ r [ h h ' + Al dxl ' { ' 

the solution can be written as 

Afi(x a ,t) = F*exp(ik a x a ), (12) 

where F\ is an amplitude of sine wave at lattice point x a and time t, k a is the wave number. 
From the Eq.tfTTl) and Eq.( fl2l) we can get F- +At = GijFj. Coefficient matrix Gij describes 
the growth rate of amplitude F\ in each time step At. The von Neumann stability condition 
is max|o;| < 1, where u denotes the eigenvalue of coefficient matrix. Coefficient matrix 
of NND scheme can be expressed as follows, 

/ At V ia At\ x Atdf? ( e ik a A Xa _ 2 + e -ik a A Xa ) 

G « = ~ T ~ -Ax^V S « + ~M 1 i^J 2 ^ 



(1 - 5) (1 - e~ ikaAxa ) , if v ia > 0; 
(l - ^) (e ikaAx * - 1) , if v ia < 0. 



1.0010- 
1.0005- 
1.0000- 

s 

£ 0.9995- 

3 

0.9990 - 
0.9985 - 
0.9980 - 



+++++ 



+ - + ;° oo % 
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— A — 


second order upwind 


1 — 
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— v— 


NND scheme 



3 

kdx 



FIG. 2: Stability analysis of several finite difference schemes. 

dfj dp <).(, or <).(, du a oj) ■ ° 2' 

There are some numerical results of von Neumann stability analysis by Mathematica. 
Abscissa is kdx, and ordinate is Iculm^ that is the biggest eigenvalue of coefficient matrix 
Gij. 

Figure 2 shows the stability analysis of several finite difference schemes. The macroscopic 
variables are set as (p, u 1: u 2 , u 3 , T) = (1.0, 4.0, 0.0, 0.0, 1.0), the other model parameters are: 
( Cl ,c 2 ,Vo) = (4.0,12.0,4.0), dx = dy = dz = 4 x 10~ 3 , dt = r = 10~ 5 , b = 5, a = p = 0. 
In this test, the NND scheme shows better stability than the others. Figure 3 shows the 
effect of dissipation term. The variables are set as (p, u±, w 2 , w 3 , T) = (1.0, 20.0, 0.0, 0.0, 1.0), 
(ci,C2,?7o) = (20.0,60.0,20.0), and the others are consistent with the Figure 2. In the two 
cases of Figure 3, operation with dissipation term is more stable (max \u\ < 1). 

Figure 4 shows the influence of parameters c±, c 2 , i] on the stability in the absence of 
dissipation term. The macroscopic variables and the other model parameters are consistent 
with those of Figure 2. Figure 5 shows the stability effect of the three parameters, when 
there is a dissipation term. The macroscopic variables and the other model parameters are 
consistent with those of Figure 3. In Figure 4 constants Ci, c 2 and f] affect the stability 
heavily. In Figure 5 the LB is stable for all tested values of c 2 and i] . Based on these tests, 
we suggest that C\ can be set a value close to the maximum of flow velocity, c 2 can be chosen 
about 3 times of the value of c±, and r] can be set to be about 1 ~ 2 times of the value of 
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i 

< 2- 



(13) 



7 



1.02 n 



□□□□□□□ 



r, X„D U 



— □— NND scheme + without dissipation term 
— o— NND scheme + with dissipation term 



1 2 3 4 5 6 

kdx 



FIG. 3: Effect of dissipation term on numerical stability. 




FIG. 4: Influence of parameters c±, C2, f]o on stability in the absence of artificial viscosity. 
IV. NUMERICAL SIMULATION AND ANALYSIS 

In this section we study the following questions using the proposed LB model: one- 
dimensional Riemann problems, and reaction of shock wave on a droplet or bubble. 
(1) One-dimensional Riemann problems 

Here, we study two one-dimensional Riemann problems, including the problem with Lax 
shock tube and a newly designed shock tube problem with high Mach number. Subscripts 
"L" and "R" indicate the left and right macroscopic variables of discontinuity. 

(a) Lax shock tube problem 



8 




FIG. 5: Effect of c±, c 2 , rjo under the condition with artificial viscosity term. 

The initial condition of the problem can be defined: 

(p, Ul ,u 2 ,u 3 ,T)\ L = (0.445,0.698,0.0,0.0,7.928), 
(p, Ul ,u 2 ,u 3 ,T)\ R = (0.5, 0.0, 0.0, 0.0, 1.142). 

Figure 6 shows the comparison of the NND scheme and the second-order upwind scheme 
without the dissipation term at t — 0.1. Circles are for the NND scheme simulation results, 
squares correspond with the second-order upwind scheme, and solid lines are for exact so- 
lutions. The parameters are (ci,c 2 ,r]o) = (2.0,6.0,2.0), 7 = 1.4, dx = dy = dz = 0.003, 
dt = t = 10 -5 . Compared with the simulation results of second-order upwind scheme, the 
oscillations at the discontinuity are weaker in the NND simulation, 
(b) High Mach number shock tube problem 

In order to test the Mach number of the new model, we construct a new shock tube 
problem with high Mach number, and the initial condition is 

(p,u 1 ,u 2 ,u 3 ,T)\ L = (100.0,10.0,0.0,0.0,0.714286), 

(15) 

(p,u u u 2 ,u 3 ,T)\ R = (150.0,0.0,0.0,0.0,50.0). 

Figure 7 shows a comparison of the numerical results and exact solutions at t — 0.25, 
where (c 1 ,c 2 ,r]o) = (8.0,24.0,8.0), 7 = 1.4, dx = dy = dz = 0.01, dt = r = 10~ 5 . The 
Mach number of the left side is 10 (Ma = u/^jT = 10/y/lA x 0.714286), and the right is 
(Ma = -u/VtT = 0). Successful simulation of this test shows the proposed model is still 
likely to have a high stability when the Mach number is large enough. 

9 




FIG. 6: Numerical results and exact solutions for Lax shock tube at t = 0.1. 
(2) Reaction of shock wave on 3D bubble problem 

The proposed model is used to simulate interaction of a planar shock wave with a bubble 
or droplet. The shock wave is moving from the right to the left. Initial conditions are (a) 



I 



(1,0,0,0,1), pre — shock, 

(p,ui,u 2 ,u 3 ,p) \ Xt y t0 = { (2.66667,-1.47902,0,0,3.94406), post - shock, (16) 

(0.1358,0,0,0,1), bubble, 



and (b) 



(p,u 1 ,u 2 ,u 3 ,p) \ x>vfl = 



(1,0,0,0,1), pre — shock, 

(2.66667,-1.47902,0,0,3.94406), post - shock, (17) 
(4.1538,0,0,0,1), bubble. 



The corresponding shock wave Mach number is 2.0, (Ma = (D — u) / V7T = (2.36643 — 
0)/Vl-4 x 1, where D = 2.36643 is the wavefront velocity). 

The domain of computation is (0 : 301,0 : 81,0 : 81). Initially, the bubble or droplet is 
at the position (200,40,40). In the simulations, the right side adopts the values of the initial 
post-shock flow, the extrapolation technique is applied at the left boundary, and reflection 
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FIG. 7: The numerical and exact solutions for high Mach number shock tube at t = 0.25. 



conditions are imposed on the other four surfaces. Specifically, at the right side, 

p(NX + 1, iy, iz) = p{NX, iy, iz) = 2.66667, 
T(NX + 1, iy, iz) = T(NX, iy, iz) = 1.6875, 
u^NX + l,iy,iz) = m(NX,iy,iz) = -1.47902, 

u 2 (NX + 1, iy, iz) = u 2 (NX, iy, iz) = 0, 

Us(NX + 1, iy, iz) = Us(NX, iy, iz) = 0, 

where ix (or iy, iz) is the index of lattice node in the x- (or y-, z-) direction, and ix = 0, 
1, • • • , NX + 1 ( iy = 0, 1, • • • , NY + 1; iz = 0, 1, • • • , NZ + 1). At the left side 
p{l,iy,iz) = 2p(2,iy,iz) — p(3,iy,iz), p(0,iy,iz) = 2p(l,iy,iz) — p(2,iy,iz), temperature 
and velocity components have the same form. Finally we take the upper surface as an 
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example to describe the reflection conditions. 



p{ix,NY + l,iz) 
T(ix, NY + l,iz) 
m(ix,NY + l,iz) 
U2{ix, NY + 1, iz) 
u^{ix, NY + 1, iz) 



= p{ix, NY — 1, iz), 

= T(ix, NY - l,iz), 

= Ui(ix, NY — 1, iz), 

= — u<i{ix, NY — 1, iz), 

= u^iix, NY — l,iz), 



p(ix, NY, iz) 
T(ix, NY, iz) 
Ui{ix, NY, iz) 
u 2 (ix, NY, iz) 
Us(ix, NY, iz) 



p{ix, NY — 1, iz), 
T(ix, NY — 1, iz), 
Ui{ix, NY — 1, iz), 
0. 

Us(ix, NY — l,iz). 



Parameters are as follows: (ci,c 2 ,?7o) — (2.0,6.0,4.0), 7 = 1.4, dx = dy = dz = 0.001, 
dt = t = 10~ 5 . Figure 8 and Figure 9 show the density iso-surfaces of bubble or droplet, 
where Figure 8 is for the process with initial condition (Tl6|) . and Figure 9 is for condition 
f fT7|) . Figure 10 shows the density contours on section z = 40, where (a) and (b) correspond 
to the processes of Figure 8 and Figure 9, respectively. The simulation results are accordant 
with those by other numerical methods 24j, |25[ and experiment 26 1. 



V. CONCLUSION 



We proposed a highly efficient 3D LB model for high-speed compressible flows. The 
convection term in Boltzmann equation is solved with the finite difference NND method, 
additional dissipation term is introduced to match the more realistic kinetic viscosity and 
to be more stable in numerical simulations. Model parameters are controlled in such a way 
that the von Neumann stability criterion is satisfied. The model can be used to simulate 
flows from subsonic to supersonic flows, especially supersonic flows with shock waves. 
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y 

FIG. 8: Density iso-surfaces of a low density bubble at t = 0.0, 0.1, respectively. 
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